선형 회귀 비용 함수의 구현

2017/02/08

n차 방정식의 미분을 취하면 다음과 같이 구현할 수 있다. numpy, pandas 쪼렙이라 순수 파이썬과 뒤죽박죽이다. R^2를 보고 차수 증가를 멈춘다. 비용함수가 컨벡스이므로 미분을 취해서 최소일 때 쎄타를 취하면 된다.

#!/usr/bin/python
import matplotlib
matplotlib.use("TkAgg")

import copy
import matplotlib.pyplot as plt
import math
import numpy as np
from scipy.stats import norm
import scipy.stats as stats
import random
from sklearn import preprocessing
from sklearn.metrics import r2_score

def frange(s, e, jump):
    r = np.array([], dtype='f')

    while s <= e:
        r = np.append(r, s)
        s += jump
    return r

min_max_scaler = preprocessing.MinMaxScaler()

orgx = frange(-1, 1, 0.01)

xx = np.array([], dtype='f')
yy = np.array([], dtype='f')
orgy = np.array([], dtype='f')

for x in orgx:
    e = random.random() * 0.1  

    if random.random() < 0.1:
       e *= 2

    if random.random() > 0.5:
        xx = np.append(xx, x + e)
    else:
        xx = np.append(xx, x - e)

yy = np.array([], dtype='f')
for x in orgx:

     y = 10 + 4 * x
#     y = 10 + 4 * x + 3 * pow(x, 2)
     y = 10 + 4 * x + 3 * pow(x, 2) + 8 * pow(x, 3)
     orgy = np.append(orgy, y)

     e = random.random() * 0.1

     if random.random() < 0.1:
         e *= 2

     if random.random() > 0.5:
        y += e
     else:
        y -= e

     yy = np.append(yy, y)

yy = min_max_scaler.fit_transform(yy)     

orgy = min_max_scaler.fit_transform(orgy)     

wi = 2
r2 = 0

while 1:

    w = [0] * wi
    minv = [100] * len(w)
    xxx = [0] * len(w)
    alpha = 5
    iter_flag = True

    while iter_flag:

        for i in range(len(w)): 
            temp = 0

    #        if minv[i] == -1:
    #            continue

            for x, y in zip(xx, yy):
                for j in range(len(w)):
                    xxx[j] = pow(x, j)

                for j in range(len(w)):
                    temp += w[j] * xxx[i] * xxx[j]

                temp -= xxx[i] * y

            w[i] = w[i] - alpha * (temp / (len(yy) * len(w)))

            if minv[i] == 100:
                minv[i] = abs(temp)

            else:
                if minv[i] > abs(temp):
                    minv[i] = abs(temp)
                else:
                    minv[i] = -1 

        for i in range(len(w)):
            if minv[i] != -1:
                iter_flag = True
            else:
                iter_flag = False

    yy2 = np.array([], dtype='f')

    for x in xx:
        temp = 0
        i = 0
        for v in w:
            temp += v * pow(x, i)
            i += 1

        yy2 = np.append(yy2, temp)

    rs = r2_score(yy, yy2)

    print rs

    if r2 == 0:
        r2 = rs
    else:
        if r2 > rs:
            break
        else:
           r2 = rs
           prevw = copy.copy(w)

    wi += 1 

yy2 = np.array([], dtype='f')

for x in xx:
    temp = 0
    i = 0
    for v in prevw:
        temp += v * pow(x, i)
        i += 1

    yy2 = np.append(yy2, temp)

plt.plot(orgx, orgy)
plt.scatter(xx, yy, s=1)
plt.scatter(xx, yy2, s=10)
#plt.plot(xx, yy2)
plt.show()

[t:/] is not "technology - root". dawnsea, rss